%% Input dataFolder and filenumbers containing evoked epscs PPRcc macro (Not randomized ISIs).
% dataFolder='7_18_17_05';
% firstFileNumber=226;
% stimInt=20;
% trials=20; 
% 
% path1='C:\Users\tomas\Documents\1. Harvard\2017 - Spring\Regehr lab rotation\Data\2017\';
% path2=['0' dataFolder(1:2) dataFolder(6:7) '\' dataFolder(1:7) '\' dataFolder '\'];
clear all
path=['C:\Users\Tom�s\Desktop\Candelabrum Cell paper\Candelabrum Cells\Mossy Fiber input\11\']
firstFileNumber=526;
stimInt=20;
trials=10; 


ISI=[10 20 30 50 100 500];
Fs=1e5; %100kHz sampling frequency
Dx=1/Fs;
tracesMat= NaN(160040,trials,numel(ISI));

for l= 1:numel(ISI); % ISI(1:end);
    
    for m=1:trials;
       
        IBW(l).files(m)= IBWread([path 'ch0_' num2str(firstFileNumber +(((l-1)*trials)+(m-1)))  '.ibw' ]);
        %trial(l).traces(m).D = IBW(l).files(m).y;
        tracesMat([1:numel(IBW(l).files(m).y)],m,l) = IBW(l).files(m).y; % create a 2D matrix where each column contains the average of each ISI)
    end

end

 
avgMat1= nanmean(tracesMat,2);  %average among columns (avg of trials).
squeeze(avgMat1);  %remove singleton dimension
avgMat=avgMat1([1:120000],:);  %trim matrix to make all traces the same lenght.
t= (0:(size(avgMat,1)-1))*Dx; %time vector using the first dimension of the trimmed average matrix.
baseline= mean(avgMat([round((0.2*Fs)):round((0.4*Fs))],:));
avgMatBS= bsxfun(@minus,avgMat,baseline); %Baseline substraction for each avg ISI.
avgMatBS(avgMatBS>3)=NaN; %remove artifact by replacing all values larger than 3pA with NaN (may need some work but for now it's fine);

avgMatBS(abs(diff(avgMatBS(:)))>2)=NaN;

% plot the average traces.
figure
plot(t(0.5*Fs:1.2*Fs),(avgMatBS([0.5*Fs:1.2*Fs],:)));
xlim([0.5 1.2]);
title([path(end-10) '.' path(end-8:end-7) '.' path(end-5:end-4) '  Cell ' path(end-2:end-1) ', 10 trial average, Stim intensity = ' num2str(stimInt)  ' �A']);
xlabel('Time (s)');
ylabel('I (pA)');
% saveas(gcf,[path1 path2 path2(end-10) '.' path2(end-8:end-7) '.' path2(end-5:end-4) '  Cell ' path2(end-2:end-1) ', 10 trial average, Stim intensity = ' num2str(stimInt)  ' �A''.tif']);


tracesMat=tracesMat([1:120000],:,:);
for l=1:numel(ISI);
        tracesMatBS(:,:,l)= tracesMat(:,:,l)-baseline(l); %this loop allows for removing the baseline from all the traces.
end

for l=1:numel(ISI);
        for j= 1:trials;
            epscMat(j,1,l)= min (tracesMatBS([0.6*Fs:0.605*Fs],j,l)); % size of epsc1 goes in column 1.
            epscMat(j,2,l)  = min (tracesMatBS([round((0.6+(ISI(l)/1e3))*Fs):round((0.605+(ISI(l)/1e3))*Fs)],j,l)); %size of epsc2 goes in column 2.
        end
end
            
epscMat(epscMat>-15)=NaN; % remove failures ( Need to find a way to do it using mad?)              
PPR= rdivide(epscMat(:,2,:),epscMat(:,1,:));
PPR=squeeze(PPR);


PPR_std= nanstd(PPR);
PPR_mean= nanmean(PPR);
PPR_SEM=PPR_std/sqrt(20);

figure;

plot(ISI,PPR_mean);
hold on;
errorbar(ISI,PPR_mean,PPR_SEM,'b')
xlabel('ISI (ms)');
ylabel('PPR');
title([path(end-10) '.' path(end-8:end-7) '.' path(end-5:end-4) '  Cell ' path(end-2:end-1) ' Mean PPR +- 1 Std, ' num2str(trials) ' trials']);
hold on;
plot([0 600],[1,1], 'k');
% saveas(gcf,[path1 path2 path2(end-10) '.' path2(end-8:end-7) '.' path2(end-5:end-4) '  Cell ' path2(end-2:end-1) ', PPR, Stim intensity = ' num2str(stimInt)  ' �A''.tif']);


avgEPSC= NaN(6,2);
 for l= 1:numel(ISI);
            avgEPSC(l,1)= min (avgMatBS([0.6*Fs:0.605*Fs],l)); % size of epsc1 goes in column 1.
            avgEPSC(l,2)= min (avgMatBS([round((0.6+(ISI(l)/1e3))*Fs):round((0.605+(ISI(l)/1e3))*Fs)],l)); %size of epsc2 goes in column 2.
 end
        
avgPPR= rdivide(avgEPSC(:,2),avgEPSC(:,1));
figure;
plot(ISI,avgPPR);


normAvgMat= bsxfun(@rdivide,avgMatBS,avgEPSC(1:6));
figure;
normAvgMat2=-1*normAvgMat
plot(t(0.5*Fs:1.2*Fs),(normAvgMat2([0.5*Fs:1.2*Fs],:)));
xlim([0.5 1.2])
title('normalized avg epscs')
ylabel('PPR');





%% load NBQX wash in and plot it
ntrials=100;
firstWashin=526
for l= 1:ntrials; % ISI(1:end);
    
    
       
        IBW2(l).data= IBWread([path 'ch0_' num2str(firstWashin +l-1)  '.ibw' ]);
       IBW2Mat(:,l)=IBW2(l).data.y-mean(IBW2(l).data.y(1:1000));
end

%remove stim artifact
a=(abs(diff(mean(IBW2Mat(:,[80:100]),2))))>2;

IBW2Mat(a,:)=NaN;
time=[1:size(IBW2Mat,1)]*IBW2(1).data.dx*1000; % in ms
figure;plot(time,mean(IBW2Mat(:,[1:20]),2),'k');
hold on
plot(time,mean(IBW2Mat(:,[80:100]),2),'r');
axis off

obj=scalebar('XLen',1,'XUnit','ms','YLen',20,'YUnit','pA')

figure;plot(IBW2Mat(mean(50220-5:50220+5),:))


